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We study the Hawking radiation in field theories which break Lorentz invariance via dissipative 
effects above a certain energy scale. We assume that the additional degrees of freedom which 
cause dissipation are Gaussian and freely falling. The asymptotic spectrum and the correlations 
are extracted from the anticommutator of the radiation field. The singular behavior of the Green 
function found for relativistic fields as one point crosses the horizon is completely suppressed by 
dissipation. Yet, when the dissipative frequency scale is much larger than the surface gravity of the 
black hole, we show that the asymptotic observables acquire their standard (relativistic) vacuum 
expectation values. We explicitly compute the effects of dissipation on the spectrum and on the 
nonseparable character of the correlations when varying the dissipative scale, the extension of the 
near-horizon geometry, and the temperature of the environment. 
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I. INTRODUCTION 

Analogue black hole radiation [1, 2] is on the verge of experimental verification, with an abundance of media now 
being studied [3]. These analogue systems mimic Lorentz invariance over large distances, while providing a natural 
mechanism for the breaking of Lorentz invariance on microscopic scales. They can thus also be used to model and test 
possible modifications of phenomena sensitive to ultrahigh-energy physics and quantum gravity effects, see e.g. [4- 
6]. High-frequency dispersion - that is, a frequency-dependent modification of the wave group velocity - has been 
well studied theoretically [7-10]. On the one hand, the robustness of Hawking radiation has been confirmed in this 
scenario when the dispersive frequency scale A is much larger than the surface gravity k determining the black hole 
temperature. On the other hand, the properties of the spectral deviations due to high-frequency dispersion have now 
been computed. We can thus expect that forthcoming experiments will observe the modified spectrum rather than 
the relativistic one computed by Hawking. In fact, the recent observations made in a water tank [11, 12] and in an 
atomic BEC [13] can only be explained by taking dispersive effects into account [14-16]. 

As yet, however, the related phenomenon of dissipation - the loss of energy to an environment - is much less 
understood, even though many systems, such as polaritons [17, 18], water waves [11, 12] and sound waves [19, 20], 
all display dissipation. 1 It is thus imperative to study theoretically the impact of dissipation on power spectra and 
correlations. It is not clear that we can infer robustness under dissipation from robustness under dispersion, since 
dissipation alters the past propagation of the outgoing modes even more drastically, and we might therefore expect that 
the deviations due to dissipation could differ considerably from those due to dispersion. Indeed, in a recent theoretical 
study of a modulated Dynamical Casimir Effect [23] , it was found that a relatively small rate of dissipation can have 
a significant effect on the properties of the emitted radiation, in particular on the entanglement of the pairs produced, 
i.e., on the nonseparable character of their state, which is the property that allows us to distinguish the stimulated 
(classical) effect from the spontaneous effect due to vacuum fluctuations [24-26]. 

In this paper, our first aim is to provide a general formalism which could be adapted so as to properly describe 
the required specific aspects of any particular system. To preserve unitarity of the system as a whole, dissipation 
is introduced through coupling to an environment. Its effects are twofold: energy is lost to the environment, but 
also gained from the environment through stimulation. Using a simple model for the environment,which allows us to 
choose how the dissipative rate varies with momentum, we shall analyze the propagation of a dissipative field in an 
analogue black hole flow. 

Our second aim is to study the spectrum of the emitted radiation and the strength of the long-distance correlations 
in order to identify the relevant parameters which govern the deviations due to dissipation. As expected, the UV 
dissipative scale Adi ss and the temperature of the environment are both important. More surprisingly, the parameter 
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1 It has also been proposed that quantum gravity could give rise to dissipative effects in the near vicinity of a black hole horizon [21]. We 
conjecture that the chaotic behavior recently discussed in [22] should effectively engender dissipation which will damp 2-point functions 
in the near-horizon region in a manner similar to that studied in the present work. 
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which fixes the spatial extension of the near-horizon geometry also plays a crucial role. When establishing this, we 
shall flesh out a prediction of [27] , that the effects caused by dissipation in the UV sector should be similar to those 
in de Sitter space since the near-horizon region can be mapped onto a portion of de Sitter space, with the surface 
gravity k identified with the Hubble constant H. 

The layout of this paper is as follows. In Section II we describe the general settings of the system we are considering, 
through the action and the resulting equations of motion. In Section III, we construct the retarded Green function 
for the radiation field in a stationary inhomogeneous situation. In Section IV, we calculate the anticommutator which 
encodes the state of the field, i.e., the spectrum and correlations between emitted quasiparticles. We then conclude 
with Section V. In appendices, we detail the calculations presented in the main body of the text. In particular, in 
Appendix B, we give the general structure of the stationary retarded Green function when the effective dispersion 
relation is a polynomial in the wave vector. 


II. SETTINGS 

In this section, we describe the structure of the models we are considering, as well as the assumptions and simplifi¬ 
cations we make for ease of calculation. We shall follow the procedure and philosophy outlined in [28]. In particular, 
the specifics of the environment and its coupling to the field do not matter so much as their desired effects once the 
environment degrees of freedom have been traced over. We therefore choose the simplest possible system compatible 
with our desire for a quantum field which exhibits dissipative effects at large momenta, and which thus reproduces 
approximate Lorentz invariance in the infrared regime. We shall mainly work with a dissipative rate which grows like 
the square of the spatial momentum, T ~ P 2 / A^iss- The dissipative UV scale Adi ss will be treated as an adjustable 
parameter to see how increasing dissipation, i.e., reducing Adi SSJ affects the observables. 


A. Metric and dispersion relation 

We work with a 1 + 1-dimensional system, and assume that the field propagates in a stationary effective metric 
which takes the form 


ds 2 = — c 2 dt 2 + {dx — v dt ) 2 . (1) 

The analogy with condensed matter systems is through the interpretation of v{x) as the local flow velocity of a moving 
medium and c as the low-frequency group velocity of the waves with respect to the medium [1, 3]. Note that K = dt 
is a Killing vector field. Its squared norm is equal to K 2 = g tt = —(c 2 — u 2 ), so there is a Killing horizon where |u| = c. 
We shall assume v < 0, describing a flow to the left, and that the magnitude of v increases in the direction of flow; 
therefore, the point where v = —c corresponds to a black hole horizon. We shall also assume c is constant, and we set 
it to 1. Note that Eq. (1) is a generalization of the (radial part of the) Schwarzschild metric in Painleve-Gullstrand 
coordinates, where v oc —1/^/x in accordance with the Einstein field equations; here, v is left arbitrary, to be treated 
as a fixed background. 

Dissipative and/or dispersive properties pick out a preferred frame [29] - typically, that in which the medium is 
at rest - which does not coincide with the stationary frame ( t , x ) since the flow velocity v ^ 0. Geometrically, the 
preferred frame is associated with a unit timelike vector field , and in 1 + 1 dimensions this also fixes (up to a sign) 
the unit spacelike vector field orthogonal to wT. If these are indeed to be associated with the rest frame of the 
medium whose flow velocity is v, then we have explicitly 

u = d t + v d x , s = d x . (2) 

These vectors are shown in Figure 1. For a given momentum p „ = (— to,k), they define the “proper” frequency and 
momentum 


S2 = —u^p^ = oj — v k , P = s^p^ = k. (3) 

The proper frequency D is thus related to the Killing frequency w simply by a Doppler shift. 

Following [27], the high-frequency behavior of the field of interest can be described phenomenologically by the 
dispersion relation 


n 2 + 2iT{P 2 )Q = fh 2 +P 2 + f{P 2 ). 


( 4 ) 
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Equation (4) gives the effective (“dressed”) dispersion relation of the held, once the environment degrees of freedom 
have been traced over. The functions T(P 2 ) and /(P 2 ) describe the deviations from the Lorentz invariant case, 
corresponding, respectively, to dissipative and dispersive effects. For simplicity, we shall set the effective mass and 
dispersion to zero, i.e. m 2 = 0 and /(P 2 ) = 0. That leaves only T(P 2 ), the dissipative rate measured in the preferred 
frame. We shall principally study the case in which T(P 2 ) = 7 P 2 for some constant 7 > 0, for this accords with 
approximate Lorentz invariance when f l and P are sufficiently small, i.e. when fl ~ P <C I /7 where 1/7 defines the 
dissipative UV scale Adi ss - However, another case of interest is that with P-independent T, relevant to systems such 
as polaritons [17]; this is briefly studied in Appendix F. 


B. Action 


The system as a whole contains two subsystems: the field of interest, <j>, and the environment field, T. For simplicity, 
we assume these to be scalar fields. The action of the system takes the form 

S = + P\]> + Pint ■ (5) 

Sj, is the action of the free <fi field. It has the form 

s 4> = \ J dtdx j[(<9 t + vd x ) cj)] 2 - [d x cl)} 2 - m 2 (j) 2 - $ f (-S 2 ) </>} . ( 6 ) 

This is just the standard 1 + 1-dimensional relativistic action for a scalar field in the metric (1), except for the last 
term in /(— d 2 ) which accounts for dispersion. The mass m and dispersion /(P 2 ) are “bare”, i.e. what they would 
be in the absence of interactions, and as such appear here without tildes. 

S-i ( is the action of the free environment field, and takes the form 


where we have defined 


S^, = — / dtdx / dq 


2 


dr^tq 



1 

T/ltt 


(d t + v d x ) y/\v\ = 7 (9* + d x v) + 7 (9t + vd x ) 


( 7 ) 

( 8 ) 


The environment is thus modeled as a continuous collection of “independent” (in the non-interacting case) oscillators, 
labeled by the continuous parameter q which is proportional to their frequency. (This parameter can be thought of 
as the wavevector along an extra dimension, as, for example, in the case of a radiating atom [30].) A dense set of 
environment modes is necessary to get a truly dissipative equation for (/>■ by contrast, a discrete set of environment 
modes would allow energy to oscillate back and forth between the 4> field and the environment. The kinetic term in 
(7), which symmetrises the placement of v and d x , has been chosen so as to yield a local field equation for <f> when 
the flow v is not homogeneous (see [27, 28]). It can be seen that the oscillators of the environment are freely-falling, 
at rest with respect to the moving medium (since they carry no spatial momentum P). The freely-falling trajectories 
along which they propagate are illustrated in Figure 1. 

Finally, Pint describes the interaction between <f> and the environment: 


Pint = / dtdxgd x (j) / dqdr'S’q 


( 9 ) 


Again, d T of Eq. ( 8 ) appears acting on \J/ q in order to yield a local field equation for </>. g is a coupling constant, 
which could depend on space and time. The non-negative integer n describes how the interaction depends on the 
proper momentum P of the cj> field: the dissipation rate T is found to be proportional to P 2n . In the main body of 
this paper, we shall principally consider the case n = 1; we focus briefly on the case n = 0 in Appendix F. 


C. Field equations 

The action described by Eqs. (5)-(9) leads to the following coupled field equations: 

[(9* + d x v) (d t + vd x ) - S 2 +m 2 + f (-<9 2 )] = (~d x ) n ^ g J dq d T ^ , 

d 2 + ( 7r «) 2 j ^ q = -d T {gd™<t>) . 


( 10 a) 


( 10 b) 
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FIG. 1. Environment and field characteristics in space-time: The dotted curves are the freely-falling trajectories, 
integral curves of the timelike unit vector field u; in the moving medium analogy, these are at rest with respect to the medium. 
They are equally spaced in position at time t = 0, and null cones have been shown for each. The two solid curves are outgoing 
null geodesics, which propagate to the right in the freely-falling frame. Note that these tilt over with x, becoming parallel to 
the Killing vector dt at x = 0; this is the Killing horizon. On its left, the flow is supersonic (|v| > c), and the outgoing geodesic 
is moving towards the left in the stationary frame. 


We commence their solution by writing a general solution for 4/: on multiplication by vfuj, Eq. (10b) can be written 
d 2 + ( 7 xq) 2 = — d T 1 from which it emerges that, at each x , acts as a driven harmonic 


oscillator. Its solution can be straightforwardly written as 




1 


VW\ J-c 


dr 1 


f sin [n q (r — r') 
7 r q 


-d r U\v\gd: 


( 11 ) 


where ’I'lj is the undriven part, i.e. a solution of the homogeneous form of Eq. (10b). When substituting the solution 
(11) into the right-hand side of Eq. (10a), the resulting equation is greatly simplified by the identity 


5 . 


dq 


sin [7 t q (t — t')\ 


7 r q 


= 5(t — t') . 


( 12 ) 


It is precisely this identity that motivates the appearance of d^q in Slj and Si nt of Eqs. (7) and (9), for it leads to 
a field equation for cf> which is local in time 2 . Still, the resulting equation can be further simplified by taking n = 0 
or 1, assuming the constancy of < 7 , and adjusting the bare mass m 2 and dispersion f(P 2 ). The derivation appears in 


2 Had we not made this simplifying choice, the wave equation (13) would have been an integralo-differential equation, containing a non- 
local kernel encoding the damping, and not a standard PDE. In homogeneous backgrounds, this can still be handled in Fourier space, 
but it becomes very messy in inhomogeneous spaces. 
















(13) 


Appendix A; here we simply state the final result: 




d t + y {-id x ) + d x (v + 1 ) 


dt + y {-id x ) 2n + (v - l)^ 


(~d x ) n J, 


where we have defined J = g f dqd T x 

As shown in Appendix A, there are two choices for the bare mass and dispersion yielding a factorized wave equation 
of the form (13), which differ in the ordering of v + 1 and v — 1. (The particular ordering in Eq. (13) is chosen for later 
convenience when g is ^-dependent; this is particularly relevant for the case n = 0 (see App. F), where it is useful to 
“switch off” the coupling asymptotically.) However, during our main focus on the n = 1 case, we shall assume that g 
is constant; then Eq. (13) becomes independent of the ordering of v + 1 and v—1. The homogeneous solutions are 
then linear combinations of solutions of the two lower-order equations 



y (- id x f n + (v ± 1 ) d x 


= 0 . 


(14) 


When g = 0, these describe waves whose group velocity in the (t,x) coordinate system is v ± 1; that is, since v is 
the flow of the medium, the waves propagate at speed 1 with respect to the medium in one of two directions. Those 
moving to the right (left) with respect to the medium will be referred to here as u- (v-) modes, respectively. Since 
we assume v < 0 , only v + 1 can vanish at some point, and hence it is principally the u-modes that experience the 
effects of the black hole horizon, which splits them into two disconnected sectors (they correspond to the outgoing 
null geodesics of Fig. 1). We shall thus focus mainly on the u-modes. It is true, however, that the presence of the 
source term on the right-hand side of Eq. (13) induces some higher-order coupling between the u- and v-modes; we 
briefly examine the strength of this coupling in Appendix E 3. 


III. DISSIPATIVE FIELD MODES 

In this section we aim to construct the retarded Green function for the </> field. This will allow us, in turn, to solve 
the sourced wave equation (13) and go on to characterize the state of the <f> field for a given state of the environment. 
The treatment here is quite general; the interested reader will find explicit expressions for the field modes and the 
Green function in a simple black hole profile in Appendix C. 


A. Retarded Green function 

The retarded Green function describes the response of the field to a delta impulse. It obeys 

□t.xGret (t, x; t', x') = S(t - t')S(x - X') , (15) 

where G\t, x is the differential operator of Eq. (13) acting on unprimed coordinates. Since this operator is f-independent, 
G Te t(t,x',t' ,x') can only depend on the time difference t — t'. It is then convenient to take the Fourier transform: 

/ +oo 

G^(i,s;^®')e < "< t - t '>d(t-i'), (16) 

-OO 

upon which Eq. (15) becomes 

n“G? et (x,x') = 6(x-x'), (17) 

where is obtained from G\t, x by substituting —iui for all occurrences of the time derivative operator d t . 

Note that, except for singular behavior along the diagonal x = x', G xet (x, x') as a function of x obeys the source-free 
form of the wave equation, and can thus be decomposed into stationary solutions of 

□TW=0 . (18) 

Asymptotically in |x| - where we shall assume v approaches constant limiting values for both x —f ±oo - G xet (x, x') 
must be a linear combination of plane waves which decrease in amplitude towards infinity. But the coefficients of 
these plane waves will generally depend on x', and so it is useful to consider how G^ et (a;, a;') behaves as a function of 
x'. Postulating that it obeys an equation of the same form as Eq. (17): 
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where the sign of ui on the operator is chosen for convenience and reflects the fact that the Green function decomposes 
into products of functions in x and x ', we necessarily have 



G re • Gret( x l,x') 


G? et (x )Xl )- □"G£ t (z ll *')) =0, 


( 20 ) 


and hence that “ is the “dual” operator defined by 




( 21 ) 


The dual operator is obtained from by flipping the sign of 7 , so \3 X U is obtained by flipping the sign of the 


entire complex frequency u — iT. Except at x = x', 
of 


G^ et (x,x') can, as a function of x', be decomposed into solutions 


□-,",rv) = o. 


( 22 ) 


Note that when 7 = 0, = □£, so that the notion of the “dual” solutions becomes redundant and the functions of 

x' in Eq. (22) are simply solutions of Eq. (18) with reversed frequency (i.e. they are complex conjugate solutions). 

Overall, G“ t ( x, x') can be written as a sum of products <jr u {xf), where the coefficient of each product depends 
only on the ordering of x and x’. Moreover, when the system is homogeneous in space (i.e. v and 7 are independent 
of x) : the entire Green function can only depend on the interval x — x'. The general case is treated in Appendix B. 
Here we only consider Eq. (13) with n = 1 and g constant. Defining 7 = g 2 / 2, the dissipative rate is T = 7 P 2 . In 
this case, the Fourier transform obeys 


-iu> + 'y{-id x ) 2 + d x (v + 1) -uo + 7 (~id x f + (v - 1) d x G“ et (x, x') = S(x - x') 


(23) 


and for a homogeneous background, we can Fourier transform also in space to get 

[—iuj + 7 k 2 + ik(v + 1 )] [—iuj + 7 k 2 + ik(v — 1)] G“ t (fc) = 1. (24) 


The roots of the first factor in (24) are solutions of the it-mode differential operator, and those of the second factor 
are solutions of the u-mode equation. We can use this property to label the roots. Hence, we can write 


G^t(A) = 


1 

7 2 (*; ~ k uR ) (k - k uL ) (k - k vR ) (k - k vL ) ' 


(25) 


The additional indices R and L indicate the sign of the imaginary part of the wave vector: it is R ( L) when > 0 
(Sfc < 0 ), so that k R () modes decay towards the right (left) side. 

By inverting the Fourier transform in /c, the R roots contribute when x > x', and vice versa for the L roots; so we 
have 


Gret.hom^a/) = i0(x — x 1 ) 


0 ik uR (x-x') 


a ik v R (x-x’) 


2(w - infk 2 uR ) 2(w - in/kl 


—i9(x' — x) 


x uR 
0 ik uL {x-x') 


vR) 

0 ikv L (x-x') 


2 (w - ijk 2 L ) 2 (w - in/k 2 L ) 


(26) 


From this exact equation, we can read out the appropriate boundary conditions that also prevail when working in an 
inhomogeneous background flow which is asymptotically constant on both sides. The key property which should be 
implemented is that G“ et ^^{x^x') is asymptotically bounded for x —> ±00 at fixed x', and for x' —> ±00 at fixed x. 
As a result, the Green function can be written in terms of 8 well-defined modes as 


Gy e t (x,x') =i0{x-x') 


Kr^Kr^') 

2 (w - ilkl R ) 


<I>vr(x)<I>vr(x') \ 

2 ( w -* 7 Kr) j 


—id( x' — x ) 


2(w - infk 2 uL ) 


2 (w - * 7 k 2 vL ) J ’ 


(27) 


where k uR and k vR in the denominator should be evaluated on the asymptotic right side, and similarly for the L 
roots on the left. The 4 solutions of Eq. (18) are easy to characterize. As an example, the dissipated mode <t>1f lR {x;) 
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is the unique solution which asymptotes on the right (for x —» oo) to the exponential e lkuRX with wave number 
k U R and a unit amplitude. Instead the characterization of the 4 dual modes, solutions of Eq. (22), is more subtle 
(a fuller derivation can be found in Appendix B). For instance <f>~ R (x') is the unique solution which asymptotically 
grows as exp(—ik U Rx') on the right side, but which decays on the left side. Hence it contains on the right side a 
uniquely determined superposition of two decaying exponentials. We could have normalized these modes, as done 
for relativistic and dispersive fields, in such a way that the denominators 2(w — i"/k 2 ) would not appear in Eq. (27). 
Since the modes have no sense by themselves due to the coupling to the environment, this does not matter here, and 
we can work with modes of unit asymptotic amplitude. Note that, in the limit 7 —> 0, the (conjugated) dual modes 
(jP(x') become standard ingoing modes, with the single asymptotic growing mode becoming the single incident wave. 

It is interesting to note that, since Eq. (13) has been fine-tuned so as to give no u-v coupling, u-modes (dual 
u-modes) are composed only of u-modes (dual u-modes) in both asymptotic regions. Therefore, Eq. (27) can be seen 
to split exactly into a sum of two terms, one containing only u-modes, the other only u-modes. We shall often focus 
our attention on the u-part of G^ et (x,x'), 


(*,*') = i6(x-x') 


•t’uR^uR^') 


2(w - *7 k 2 uR ) 


— i0(x' — x) 




2 (w - i^k 2 uL ) 


(28) 


When considering a black hole flow, there is an interesting switch amongst 4>“ R and 4>^ L - When the flow vanishes 
or is subsonic, the u^-mode is the dissipative version of the standard u-mode, whereas the additional root k U L is 
brought in by the dissipative term in k 4 in the dispersion relation. In fact, its imaginary part diverges as —(v + l)/q 
when 7 —► 0. When the How is supersonic and to the left (i.e. u + 1 < 0), it is now the uz,-mode that is the dissipative 
version of the standard it-mode, while k U R is the additional root with large imaginary part. This can be understood 
from the fact that, in supersonic flows to the left, the group velocity of u-modes in the lab frame is to the left, see 
Fig. 1. This switch can be studied by considering the roots of the u-dispersion relation, see Eq. (14), 

[—ioj + 7 k 2 + ik u ( v + 1)] = 0. (29) 

One sees that the standard root in one region varies continuously to become the additional root in the other region. 
The migration of the k U R root is shown in Figure 2 as u +1 varies from D (subsonic region) to —D. (When considering 
the u-modes, no such switch arises as u — 1 does not flip sign in a black hole flow to the left.) 


B. Adimensional units 


When focusing on u-modes, we shall see that the dissipative wave equation (14) defines an adimensional parameter 
which depends both on the coupling strength 7 = g 2 / 2 and the properties of the background flow. Recalling that 
c = 1, notice that 7 has dimensions of time, as it is related to the dissipative rate by T = 7 P 2 = P 2 / Adi ss - 
Since u-modes are sensitive only to the sum v + 1 we can write 

v(x) + 1 = Dh(^jpj , (30) 

where k and D are parameters while h varies monotonically between two asymptotic limiting values, with h( 0) = 0 
and h'(0) = 1. Then there is a horizon at x = 0, with “surface gravity” |u'(0)| = k. For a given profile, the parameter 
D allows us to vary the total difference in v between the asymptotic regions, or, equivalently, the size of the near¬ 
horizon region (NHR) where v{x) + 1«ki. In the limit D —> 00 , v(x) + 1 becomes exactly kx, which describes the 
Poincare patch of de Sitter space in the stationary picture [27]. 

We can adimensionalize the u-mode wave equation of (14) by defining dimensionless time and space variables: 


T = 


t 

to 1 



(31) 


By direct substitution, we get 


— dr + 
to 


-2 (~idxf + —h 
xi x 0 



0 = 0 . 


(32) 


It is natural to select to = 1/re, independent of the dissipative rate. However, there are two natural choices for 
the length scale Xo : there is a “geometrical” length Xq 60 " 1 adapted to the geometry described by h; and there is a 
“dissipative” length Xq 1ss adapted to the effective rate of dissipation in space: 


D 

K 


X 


diss 

0 



geom _ 
— 


(33) 







FIG. 2. Roots OF DISPERSION RELATION of Eq. (29): Here is shown, on the complex k U R- plane (and in logarithmic scale), the 
continuous variation of the wave vector k U R as v 4 - 1 varies between D and —D, i.e. as the position is varied from the subsonic 
to the supersonic region. k U R has positive real and imaginary parts, and is thus outwardly propagating to the right. In the 
subsonic region, as x — > oo, its imaginary part takes its smallest value, and this corresponds to the lower end of each curve. 
As x is varied back towards the horizon, its real part gradually increases, turning around at x = 0. As i is further varied into 
the left-hand (supersonic) region, the real part of k U R decreases again while its imaginary part continues to increase, so that 
the top end of each curve corresponds to a: -I —oo. The various styles of curve correspond to different values of A 2 = 7 k/D 2 \ 
0.05 (solid), 0.1 (dashed) and 0.2 (dotted). The three different groups correspond to different values of oj/k. The three light 
dot-dashed curves plot the roots at fixed v + 1 for varying u, with A 2 = 0.05. Note that the ( v + l) su b = D curve follows the 
expected pattern with Iin(fc) oc (Re(fc)) 2 ; the ( v + l) aup = —D curve has Im(fc) independent of Re(fc); and the v + 1 = 0 curve 
(i.e., at the horizon) has Im(fc) = Re(fc), and is the same regardless of the value of A 2 . 


Calling X = x/x^ om and \ = ^c/;Cq 1ss , Eq. (32) can therefore be written in the equivalent forms 


1 9t + A 2 (— idx) 2 + h (A) dx 


d T + {-id x ) 2 + l ^d x 


where we have defined the dimensionless parameter 


0 , 

0 , 


(34a) 

(34b) 


A 2 




2 


7 K 

D*' 


geom 


(35) 
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Re[2co G? etu (x,x')\ 


Re[2u) G% lu (x,x')] 




K X' 

D 



FIG. 3. Green function: Here is plotted the u part of the retarded Green function, G? etiU (x, x') as a function of the source 
coordinate x' for fixed image coordinate x. The various curves correspond to different values of A 2 = 7 k/D 2 -. 0 (dotted curve), 
0.1 (dashed curve) and 0.2 (solid curve). The columns correspond to different frequencies, with lj/k = 1/2-71 on the left and 1/2 
on the right. In the top row, the position coordinates have been adimensionalised by the “geometrical” length (the first of Eqs. 
(33)) and the Green function has been multiplied by 2 uj so as to see the effects of dissipation through a reduction with respect 
to a unit amplitude. In the bottom row x has been adimensionalised by the “dissipative” length (the second of Eqs. (33)). 
In this representation, the limiting case A 2 = 0 corresponds to working with a dissipative field in a de Sitter geometry. The 
small differences with the cases A 2 = 0.1 and 0.2 illustrate that the near-horizon behavior of G? et}U (x, x') in a black hole flow is 
essentially the same as that in de Sitter. The vertical dashed lines show the fixed position of the image coordinate: kx/D = 10 
(in the flat asymptotic region) in the top row, and X\Jn/'y = 2 (well inside the near-horizon region) in the bottom row. 


Equations (34) are both useful in different regimes: they show how position should be scaled in order to remove any 
significant dependence on A 2 . So when A 2 P 2 -C hP , as is most probable in the asymptotic regions, the dissipative 
term in Eq. (34a) can be treated as a small perturbation, and the field modes do not depend significantly on the 
particular value of A 2 . This assumption breaks down near the horizon where P becomes very large, but there Eq. 
(34b) becomes independent of A 2 since h(X\)/X « x when Ay is sufficiently small. 

Note also the different limits probed by Eqs. (34a) and (34b) when A 2 —> 0: Eq. (34a), treating the “geometrical” 
length as the fundamental length unit, sends the dissipative rate to zero so that we recover the standard relativistic 
wave equation; on the other hand, Eq. (34b), which treats the “dissipative” length as fundamental, effectively sends 
D to 00 as A 2 as decreased. In this limit, one is considering a dissipative field in de Sitter space [27], since v + 1 = kx 
for all x. 

In Figure 3 are plotted examples of (the u part of) the retarded Green function, for various dissipative rates, in 
both the “geometrical” coordinate X and the “dissipative” coordinate x- Note the effect of dissipation according 
to these two viewpoints. For fixed X with varying X\ the largest changes induced by dissipation occur around the 
horizon, where the infinitely many oscillations of the relativistic model are smoothed into an exponential tail; in the 
asymptotic region, dissipation causes a gradual decrease in the amplitude, but we can see from the position of the 
node that the phase of the function is hardly affected. On the other hand, when \ is fixed inside the near-linear 
de Sitter region and x' is varied, the behavior around the horizon possesses a well-defined (and nontrivial) limit as 
A 2 —>■ 0. Whenever A 2 1, this establishes that the behaviour in the near-horizon region is indeed very close to 
that found in de Sitter space (the deviations vanish as A 2 for A 2 —>■ 0). The plots of Fig. 3 have been obtained in a 
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particular flow where the decaying inodes can be explicitly computed. This simple case is studied in Appendix C. 


IV. STATE OF THE FIELD 

In this section we include the effects of the source term {—d x ) n J in Eq. (13) in order to calculate various observables 
characterizing the state of the <f> field. The statistical properties of the operator J are fixed by the initial state of the 
environment, which is assumed to be thermal, completely characterized by its temperature Tip. The entire information 

is thus contained in the anticommutator ^|j, j|^, which is called the noise kernel [31]. Because the field equations 

(10) are linear and because the background fields g^ and zz p are stationary, the state of the </> field will also be 
Gaussian and stationary. As a result, each w-sector can be studied separately, and its state is completely determined 
by the w-components of the anticommutator, 

G“ c {x,x') = ({(TOr),^"^^)}) , (36) 

where to is real and positive, and where 0 _ “(a:) = since the field is hermitian. 


A. Interpretation of the anticommutator 


To interpret the form of G“ c , let us first consider the w-component of the field operator in the non-dissipative model 
(using notation adopted from [26]): 

4> u {x) = a u p“ R (x) + (a.J 1 (^“(i))* + al^(x). (37) 

The amplitude operators obey the usual bosonic commutation relations, and the relativistic modes tp u are normalized 
according to the standard convention (~ e~ lut /y/2uj). For u> > 0, the ztL-mode has negative norm while the others 
have positive norm, and this is indicated by writing the zzL-mode contribution to the field operator as the hermitian 
conjugate of the same term with frequency — w, which has positive norm. This allows us to suppress the u label on 
the amplitude operators and distinguish the left and right amplitude operators by the sign in front of u>. 

By direct substitution of Eq. (37) in Eq. (36), in any stationary state, we find 


G“ c (x,x') = (2n w + 1) <p“ R (x) {<Pur( x '))* + (2n_ w + 1) {<p u £(x))* <p u £(x') + (2n” +1) <p“(x) (<p%(x'))* 

+ 2 Cu Vur{x)vZl( x ') + 2c l Vv( x )vZZ{x') + 2 d v u <p“(x) (<p“ R ( x')f + (x o x'f , (38) 


where the additional contribution (x <-> x')* only applies to the second line. We have used the bosonic commuta¬ 
tion relations that exist between the amplitude operators in order to replace the anticommutators with standard 
observables: 


riu 



Cw — 





(39) 


These observables, then, are directly related to the amplitudes of various terms appearing in the decomposition of the 
anticommutator. In turn, following [26] we can construct from them the following quantities: 


= 


n^ri-oj Co. 


A" = 




D v „ = n u 


- K 


(40) 


If the state of the system can be described classically, all of the quantities in Eqs. (40) must be strictly nonnegative. 
Quantum mechanics, however, allows and A" to be negative, bounded from below by the corresponding values of 
—n. (D^j, examining correlations between quasiparticles of the same norm, remains strictly nonnegative even in the 
quantum description [24].) States for which either or A” is negative are said to be nonseparable: the strength of 
the correlations is so strong that they cannot be described by a classical ensemble. 

We shall principally be interested in A w , which measures the “quantumness” of the pairs of w-quasiparticles which 
propagate on either side of the horizon. To extract A w , we shall use the fact that G“ c is a sum of three terms: 


G% c (x,x') = G^ u (x,x') + G^ v (x,x') + G“ 




(41) 


The first term governs the sector of interest containing the two zt-modes, the second one the sector containing the 
single z;-mode, and the last one contains the zz/u-correlations governed by c“ and In the present model, the latter 
do not vanish even though the retarded Green function is a sum of a u and a v contribution, for both sectors are 
sourced by the same environment. However the strength of the u/v correlations are very small (see Appendix E3). 
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B. Behavior of the anticommutator 


From the wave equation (13), taking the Fourier transform in time, we obtain 

n“n-,“G^(x,x') = d x d x , = d x d x ,N“(x,x'). 


(42) 


Note that the operator □ appearing in Eq. (42) is not a dual operator of Eq. (22). The solution can be built from 
the retarded Green function: 


G“ c {x,x') = J dx i J dx 2 G“ et (x,xi)G ie t(x',X 2 )d Xl d X2 N u (x 1 ,X 2 ) 

= j dx i Jdx 2 d Xl G^ et (x,x 1 )d X2 G~^(x',x 2 )N u (x u x 2 ), 


(43) 


where in the second line we have used integration by parts to move the derivatives from the noise kernel to the Green 
functions. The anticommutator is thus determined from a knowledge of the retarded Green function (calculated in 
§111 A) and of the noise kernel (calculated in Appendix D - see Eq. (D.ll)): 

J e iu{T 0 {xi)-T 0 (x 2 )) r+ °° 

N u (x 1 ,x 2 ) = -- - , / dlOqUq coth 

71 v(X\)v[X 2 ) J-oo 

p iui(T 0 (x 1 )-T 0 (x 2 )) r+oo / \ 

= — ; / dujq coth — (iv(x i) 0 Xl - iv(x 2 )d X2 ) e -«-«M* 0 -n>(*»)) ( 44 ) 

2?r ^v( Xl )v(x 2 ) 7-00 9 V2 Tv) y K X1 y 21 X2> y ; 

where in the second line we have rewritten ui q in terms of derivatives of the exponential, and where 


/ U>q \ e -iui q (T 0 (x 1 )-T 0 (x 2 )) 

V 27 vrv / 


7o(®) 



(45) 


is the proper time along a freely-falling geodesic between positions 0 (the horizon) and x. (Since v < 0, tq(x) decreases 
when increasing x.) Plugging the expression for the noise kernel back into the expression for the anticommutator, we 
can integrate aq and x 2 by parts and integrate over uj q - using the fact that the Fourier transform of the coth function 
is another coth function of inverse width, as can be shown by expanding it into a sum of poles in the complex plane 
- to obtain 


r+00 

f +00 1 

/ dx\ 

/ dx 2 < 

J —00 

7 — 00 [ 


glOJToOl) 

-v/K^i)! 


d Xl G“ et (x,xx) ■ d, 


0 —iuiTo(x 2 ) 


-v{x 2 )d X2 G-?(x' ,X 2 ) 


-d x 


piuJTQ (xi) 

\/K2qJ| 


v( Xl )d Xl G? et (x ,Xl) 


e ~iuT 0 {x 2 ) 


\ZHx 2 )\ 


V\ V ( X 2)\ 

dx 2 G~ e “{x',x 2 ) > x 7rT^coth(7rT^(ro(xi) - r 0 (x 2 ))) . (46) 


This is a useful expression in that the singularity at aq = x 2 is of the form of the simple pole ( X\ — aq) -1 , an d can thus 
be integrated by taking the Cauchy principal value. An effective way of doing this is to switch the dummy variables 
X \ and £ 2 , and take the mean of the two integrands; since coth(z) is an odd function, this has the effect of removing 
the singularity altogether, so that numerical integration can proceed. 

In Figure 4 we illustrate the integrand of Eq. (46) having re-expressed the integral using the adimensional lengths 
Xi and X 2 (see Eqs. (33)) as dummy variables. Using these variables, the limit 7 —>• 0 is finite. The plots show 
the integrand along the diagonal xi = X 21 having sent x and x 1 into the right-hand asymptotic region and extracted 
the product of two it-R-modes ^ i j(x)0“^(£ / )/(2w) 2 . We clearly see that x is the appropriate integration variable, 
with variations in A 2 and ui/k only slightly deforming the shape of the integrand in these units. Note that, as well 
as the modes in x and x', two factors of l/2w have been extracted so that, as can be seen in the left plot of Fig. 4, 
the remaining integrand approaches a well-defined limit as ui —> 0. One of these factors “normalizes” the asymptotic 
plane waves, so that the coefficient of q i >ur( x ) c I } ur( x ')/2u corresponds to 2n w + 1. The second factor of l/2w shows 
that, as w -> 0, + 1/2 ~ 1/w, which corresponds to the infrared divergence of the Planck spectrum. On the right 

plot, we can verify that the integrand possesses a well-defined and non-trivial limit when A 2 —>■ 0, guaranteeing that 
the anticommutator is a well-defined function of x and x' in the relativistic limit. 

In Figure 5 are plotted examples of the u part of the anticommutator G“ c u (x,x r ) as a function of the source 
coordinate x\ x being fixed at kx/D = 10, far in the subsonic asymptotic region; the geometry is fixed at D = 0.5, 
and the environment temperature is very low: T i>/k = 10~ 2 . (This case is very close to the zero temperature limit 
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Int■ (xi = Xi) 



Int. (xi = Xi) 



FIG. 4. Integrand for anticommutator: Here is plotted the integrand of Eq. (46) re-expressed in adimensional “dissipative” 
length units \ = X\J re/7. We plot the integrand along the diagonal = X 2 , having sent x and x' into the right-hand 
asymptotic region and factored out the product of uR -modes (t>^ R {x)(j>~ R {x') / (2uj) 2 . In the left plot, A 2 = 0.1 while the various 
curves correspond to different values of ca/re: 1/207T (solid), l/2n (dashed), 1/3 (dot-dashed) and 1/2 (dotted). In the right 
plot, lo/k = 1/2-7T, while the various curves correspond to different values of A 2 : 0.1 (solid), 0.2 (dashed) and 0.4 (dotted). 
Throughout, the geometry is fixed at D = 0.5. We clearly see that the integrand is peaked in the near-horizon region with a 
half-width of about one dissipative length zj) (see Eqs. (33)). This remains true in the limit A 2 -A- 0. 


since T^/Th ~ 1/15, where Tfj = k/2t: is the Hawking temperature.) More precisely, we represent the first term of 
Eq. (41) (found by replacing G“ t with its u part G“ t u in Eq. (46)) multiplied by 2 w so as to be able to read off 
the values of the power spectrum (= 2 n u + 1) for x' = x, and of the strength of the correlation 2c u for x' = —x, 
see Eq. (38). Various values of A 2 have been used, including A 2 = 0. In this relativistic case, the state of the field 
has been taken to be the standard Unruh vacuum [32]. The two columns correspond to different frequencies, with 
w/re = 1/27r in the left column and w/re = 1/2 in the right column. 

Three important observations should be made. Firstly, on the right side of the two upper plots, by comparing 
the two curves (A 2 = 0.05 and 0.1) with the relativistic case, we can see with the naked eye that dissipation does 
not have a significant effect on the mean number of quasiparticles emitted to the right. This means that the low 
temperature fluctuations of the environment feed the radiation field (in a vicinity of the horizon approximately given 
by one dissipative length Xq 1ss of Eq. (33)) in such a manner as to bring its state to the Unruh vacuum in the 
limit A 2 —A 0. The smallness of the deviations for A 2 C 1 confirms that the asymptotic spectrum is robust against 
introducing dissipation in the ultra-violet sector. (In appendix E2, we shall examine to what extent this remains 
true when increasing the temperature of the environment above the Hawking temperature.) On the right upper plot, 
notice that the asymptotic value ~ 1 means that we are close to the vacuum for this value of the Killing frequency 
w = re/ 2, i.e., about three times the Hawking temperature Tj-j = k/2it. 

Secondly, as for the retarded Green function, see Fig. 3, dissipation smooths out the logarithmic phase singularity 
(in \x'\ w / K as |cc'| —>• 0 ) that occurs at the horizon in the relativistic case; this is most clearly seen in the second row, 
which zooms in on the near-horizon region. In this sense, the “trans-Planckian” behavior in the vicinity of the horizon 
is completely erased in our dissipative model with T = 7 P 2 . 

Thirdly, and most interestingly, unlike for the retarded Green function and despite the large deviations caused by 
dissipation in the near-horizon region, it is seen that on the other side of the horizon the anticommutator redevelops 
so as to smoothly join to the oscillatory behavior that occurs in the absence of dissipation (and dispersion). While 
the amplitude is slightly decreased with increasing dissipation, the nodes are seen to occur at the same places. As 
a result, the relative phase between the two asymptotic regions, which is encoded in the phase of the c u coefficient 
of Eq. (39), is preserved. In this we generalize the observation that n w and c u are both robust when introducing 
high-frequency dispersion [33]. In other words, the asymptotic state of the u- modes is robust against UV dissipation, 
i.e., the differences of and c u with respect to the relativistic values linearly vanish with A 2 when A 2 1. 

Yet, in order to determine the domain where the state is nonseparable, which is that where of Eq. (40) is 
negative, we need to extract the occupation numbers and the correlation term with precision. 
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Re[2a> G“ V)] Re[2w G^ cu (x,x ')] 





Re[2cj G^ cu (x,x')] 



FIG. 5. Anticommutator: We plot here the real part of the M-mode contribution to the anticommutator, G% CtU (x, x'), as 
a function of x', with x fixed at kx/D = 10, the geometry fixed at D = 0.5 and the environment temperature fixed at 
Ty/n = 1CF 2 . The dotted curves plot the non-dissipative limit (note that G^u 4 ' has a phase singularity as x’ —> 0); the dashed 
and solid curves have A 2 = 0.05 and 0.1, respectively. The two columns correspond to different frequencies: uo/k = 1/2-r (left 
column) and 1/2 (right column). The second row shows the same information as the first row, but zoomed in to the near-horizon 
region, so that effects of dissipation can be better observed. As explained in the text, up to (small) damping factors discussed 
in the next section, on the asymptotic right side, the value of (2lj) x G ac ,u gives 2 x 1 ^ + 1 whereas when evaluated on opposite 
sides it gives the real part of 2c aj e* 9 ^’ :E > where the phase 6(x,x') accounts for the propagation from x to x 1 . We clearly see 
that dissipation in the UV removes the singular behavior as x’ —> 0 of the relativistic case while preserving the asymptotic 
properties of G ac ,u on both sides of the horizon. 


C. Spectrum and entanglement for low environment temperature 


Because of the residual coupling to the environment, the modes are still slightly damped in the asymptotic regions. 
Hence, the extraction from the anticommutator of the values of n u and must either be done at a finite distance 
from the horizon, or one should explicitly switch off the interactions between the radiation field and its environment. 
As a result, the identification of n u and c u is inevitably slightly imprecise. The importance of the errors is governed 
by r as /w, where T as is the asymptotic residual decay rate. For low values of A 2 of Eq. (35), this is not a serious 
problem in our model because F as /w is equal to X 2 lo/k. Similar errors have been discussed in homogeneous time 
dependent settings in Sec. IV.B.2 of [34]. To minimize the total absorption during outwards propagation from the 
black hole, while also having a locally flat geometry, we work at k|:e| = 5 D. 

In Figure 6 we plot the spectral deviations on the subsonic side when the environment is again at a low value of 
T^/n = 10~ 2 , the velocity profile fixed at D = 0.5 and the positions of the hypothetical particle detectors fixed at 
k\x\/D = 5. The limit A 2 —> 0 thus corresponds to the relativistic limit. (We do not represent the spectrum of 
negative energy partners in the left (supersonic) asymptotic region because it behaves very similarly.) From the left 
plot, we see that for low frequencies, the relative difference in the quasiparticle numbers from their relativistic values 
tends to a constant, which translates into an effective change (specifically a reduction, since decreases) in the black 
hole temperature. At higher frequencies, we notice that the relative difference is greater, meaning that the spectrum 
falls off more rapidly than it would for a purely Planckian spectrum. 

To see to what extent this reduction is due to damping of the modes as they propagate from the horizon to the 










14 



FIG. 6. Power spectrum at low temperature: Shown here is the relative change in the number of quasiparticles from 
the Hawking thermal prediction, in the right (subsonic) asymptotic region. The temperature is fixed at Ty/n = 10~ 2 , with 
measurements made at a fixed distance from the horizon of k \x\ /D = 5. The geometry itself is fixed at D = 0.5. In the left 
column, the variable is lo/k, and the various curves correspond to different values of A 2 : 0.1 (solid curve), 0.2 (dashed curve), 
0.3 (dot-dashed curve) and 0.4 (dotted curve). The dots show the data points of the solid curve corrected for the residual 
damping outside of the near-liorizon region, which is most significant at high frequencies. In the right column, the variable is 
instead the adimensional parameter A 2 of Eq. (35), while the various curves correspond to different values of w//c: 1/20-7T (solid 
curve), 1/2-7T (dashed curve), 1/3 (dot-dashed curve) and 1/2 (dotted curve). On both plots, we clearly see the linearity of the 
deviations in A 2 when A 2 < 1. 


detector, we multiplied n^ by the correction factor exp(2 JZ,„ dx'Qk UJ (x')) which accounts for the damping from x em ‘ 
to x det '. If the final value a; det ’ is clear, and given here by nx/D = 5, some rule should be adopted for determining 
the locus of emission x em ~. We adopt the location where the diagonal of the integrand of Eq. (46) is equal to 1/5 
its maximum value, see Figure 4. The result is represented by the dots in Figure 6. From this we conclude that the 
higher damping observed at high frequency is largely due to the mode propagation, and not intrinsic to the black hole 
emission. 

From the right plot, we see that, for a fixed frequency to, the relative difference in occupation number linearly 
increases with the parameter A 2 when it is smaller than 0.1. For uj/k ~ 1, the coefficient of proportionality is of order 
1. The fact that the relative difference in due to dissipation is proportional to A 2 was anticipated in [27] on the 
basis that, firstly, the near-horizon black hole physics is similar to that occurring in de Sitter space, and secondly, 
that the relative difference in n u due to dissipation in de Sitter was found to be extremely small. (The subdominant 
deviation of the spectrum due to D at fixed A 2 is studied in appendix El.) It is therefore the departure from the 
de Sitter physics that fixes the leading spectral deviations of the black hole spectrum. It should be noticed that 
similar observations apply to dispersion effects [10], and also to the local description of the correlations (encoded in 
the anticommutator): as long as one probes the G ac in the de Sitter-like near-horizon region, there is no departure 
from the vacuum structure (both for relativistic fields, and for weakly dispersive/dissipative fields). The typical long 
distance behavior [35, 36] is gradually obtained as one leaves the de Sitter-like region [37]. 

In Figure 7 we plot the parameter of Eq. (40), the sign of which determines if the final state is entangled or 
not. We observe a strong robustness of the nonseparability of the emitted quasiparticles: no loss of nonseparability 
is observed for the values considered, though it appears that it would be lost when we consider higher frequencies, 
to /k 1, which are subject to dissipation outside the near-horizon region. The strong robustness can be understood 
as a new manifestation of the fact that the mode mixing responsible for the Hawking effect effectively occurs at low 
momentum P ~ k. In effect, strong dissipation occurs at high momenta P, before vacuum configurations are converted 
into on shell particles, as discussed above. Therefore dissipation does not significantly reduce the entanglement. When 
studying dissipative effects in cosmological settings, the importance of the time ordering of dissipation and pair creation 
was clearly established [38]. As a result, when working in de Sitter space [27], for a given field with dissipative effects 
in the deep UV sector, it was found that nonseparability can be lost for the cosmological pairs with opposite momenta 
(k, — k ) (because they are produced when dissipation is still important), whereas it can be preserved for the -u-pairs 
(w,-w). 
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FIG. 7. Nonseparability at low temperature: Here is plotted the rescaled nonseparability parameter, A y /n u of Eq. (40), 
which is bounded from below by —1; to make the differences more visible, 1 + A^/n^ has been plotted in logarithmic scale 
in the left plot. As in Fig. 6, we work with T^/k = 10~ 2 , and k\x\ /D = 5, where D = 0.5. In the left column, the variable 
is lo/k, and the various curves correspond to A 2 : 0.1 (solid curve), 0.2 (dashed curve), 0.3 (dot-dashed curve) and 0.4 (dotted 
curve). In the right column, the variable is instead A 2 of Eq. (35), while the various curves correspond to a j/k: 1/207T (solid 
curve), l/2n (dashed curve), 1/3 (dot-dashed curve) and 1/2 (dotted curve). Even though the departures from the relativistc 
case A 2 = 0 are clear, quite remarkably, in all represented cases, the two-mode states of u-quanta remain nonseparable. 


V. CONCLUSIONS 


We have studied the effects of dissipation on the properties of the quasiparticle pairs emitted by an analogue black 
hole. Since dissipation implies the absence of well-defined particle modes and their corresponding annihilation and 
creation operators, we found the properties of the radiation field through examination of its anticommutator, which 
for our Gaussian model completely determines the asymptotic state of the field. We worked with 1 + 1-dimensional 
stationary models. 

To preserve unitarity of the system as a whole, we began with an action that included a dense environment, which 
was later traced over to obtain a driven-dissipative equation for the radiation field. The kinematical term of the 
environmental degrees of freedom, and their coupling to the radiation field were fine-tuned so as to yield a simple, 
local dissipative equation. In particular, the decoupling of u- and c-modes (which is obtained for massless relativistic 
fields in two dimensions) was preserved for the dissipated modes in a model where the decay rate grows like the 
square of the momentum measured in a freely falling frame (T oc P 2 ). This allowed relatively easy computation of 
the retarded Green function, which was found to depend on only two dimensionless parameters: the frequency uj/k 
and A 2 = 'yn/D 2 , which combines the dissipative frequency Adi ss = 1/7, the black hole surface gravity n, and D , 
the extension of the near-horizon region which can be mapped onto a stationary portion of de Sitter space. As a 
result, the limit A 2 —> 0 can be conceived in two different ways: as giving the relativistic limit in a fixed metric, or 
as working in de Sitter space with a fixed dissipative frequency Adi ss - At fixed Killing frequency, we showed that 
the retarded Green function can be canonically expressed in a double mode basis containing, on the one hand, the 
dissipated modes, and on the other hand, dual modes obeying well-defined conditions on asymptotic spatial infinities. 
The origin of this writing, which should apply to all unitary models as indicated in Appendix B, comes from the fact 
that the retarded Green function itself is the unique solution which is asymptotically bounded. We then studied its 
near-horizon properties in the model T oc P 2 , and saw that the singular behavior of the relativistic theory in 6{x) x lu> ' K 
as x —> 0 was fully regularized by dissipation. 

The Green function was then used to calculate the anticommutator, the calculation of the double integral being 
done numerically. Focusing on the T oc P 2 model, the asymptotic spectra and long-distance correlations of the emitted 
pairs were calculated along with the local properties of the anticommutator. Two apparently contradictory results 
were found. On the one hand, when sending one point into the near-horizon region while keeping the other fixed far 
away, we recovered the radical departure from the relativistic case (already found when studying the retarded Green 
function) by which the singular behavior at the horizon is smoothed out. On the other hand, despite this radical 
departure close to the horizon, we find that when A 2 < 1, the anticommutator far from the horizon in both directions 
behaves essentially like that of a relativistic field in the Unruh vacuum. In other words, the strong interactions 
occurring at high momenta P 2 ~ kA^ss in the near-horizon region are such that the fluctuations of the environment 
bring the state of the radiation field to its regular ground state. The origin of this fact can be understood by analyzing 
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the behavior of the integrand of the double integral defining the anticommutator: one finds that it possesses a scaling 
property which guarantees that the decoupling limit A 2 —> 0, i.e., A<ji ss —> oo, is non-trivial and implements the 
fluctuation-dissipation theorem, namely that the state of the radiation field is the same as that of its environment. 
We then computed the leading spectral deviations due to dissipation. For a fixed value of lo/k, we found that the 
deviations principally (and linearly) depend on A 2 . This confirmed that the spatial extension of the de Sitter-like 
near-horizon region plays a crucial role in determining the deviations. 

In Appendix E 2 we verified that the asymptotic spectrum and the correlations remain as in the Unruh vacuum 
when the temperature of the environment is higher than the Hawking temperature, provided dissipation only occurs 
above a sufficiently high energy scale. The robustness stems from the large redshift which effectively suppresses 
the stimulated emission. Hence, when Adiss is sufficiently high, only the spontaneous contribution due to vacuum 
fluctuations is responsible for the asymptotic properties. This was explicitly verified by calculating the parameter 
which measures the degree of the nonseparability of the emitted pairs: we observed that the final state remains 
entangled over a wide range of environment temperatures and dissipative rates. 

In conclusion, even though we worked with fmed-tuned dissipative models, mainly characterized by an absence of 
dispersion in the limit A 2 —» 0 and decoupling of the u- and c-sectors at the level of the field modes, we claim that our 
findings are generic in character. Our principal findings concern, on the one hand, the mode structure of the retarded 
Green function at fixed w, and how it determines both the local and the asymptotic properties of the anticommutator; 
and on the other hand, the robustness of the spectrum and the correlation strength when A 2 < 1, and the properties 
of their leading deviations. These results should also apply to more complicated models which describe analogue black 
hole flows in water waves, in acoustics, and in polariton systems. Note however that the latter case is more sensitive 
to dissipation as the rate T is independent of the momentum; see Appendix F. 
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Appendices 


Appendix A: Derivation of the driven-dissipative wave equation 

In this appendix, we show how Eq. (13) is derived from Eqs. (10)-(12) plus some additional simplifying assumptions. 
Substituting Eq. (11) in Eq. (10a) and using Eq. (12) yields the following: 

(dt + 0 x v ) (d t + vd x ) - d 2 x + m 2 + / (~d 2 x ) + (-1 ) n dZgd r gd*\ <j> 

(d t + d x v ) (d t + vd x ) - + m 2 + f (-d 2 x ) + y -^-d n x g (d t + d x v) gd n x + [ -^ L dy (d t + vd x ) gd" 

= (-l ) n d%g [ dqd T *° q , (A.l) 


where in the second line we have used Eq. (8) to expand d T . It is convenient to combine the last two terms of the 
second line with the first term, but to do this we must first change the order of some of the differential operators. We 
require the following operator identities: 


d n x g (dt + d xV ) gd: = (dt + d x v) d n x g 2 d: + (d*gd x v - d x vd n x g) gd n x 

= (d t + d x v) d n x g 2 d: + (d n x +1 gv - d n x g'v - d x vdy) gd,, 
= 2 (dt + d x v) dyd: + 2 d x (dy - vd n x ) 7 s : - dyid: 
d n x g (dt + vd x ) gd: = d n x g 2 d n x (d t + vd x ) + d n x g (vd x gd n x - g^vdj 

= d n x g 2 d: (dt + vd x ) + dy ( vg d: +1 + vg 'd: - gd n x vd x ) 
= 2dyd: (dt + vd x ) - 2(9 "7 (dy - vd n x )d x + dyyd: 


j 

X > 


/ rp 
x > 


(A.2a) 


(A.2b) 
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where in the last lines we have defined 7 = g 2 / 2, and a prime denotes the ^-derivative. Using Eqs. (A.2) to replace 
the last two terms of Eq. (A.l) and rearranging, we find 


[(& + d x v + (—l) n dydy (d t + vd x + (-1 ) n d£ 7 d£) - 9 2 
+ m 2 + f (- a l) - (ayay 2 + (■ -I) n d x (ay - vay 7 a n x - (-i) n ay (ay - vay a x 


(~d x ) n J(t, x), (A.3) 


where J(t,x) = g f dqdy Q q . 

Equation (A.3) is quite general: v and 7 are arbitrary functions of x , n is an arbitrary (nonnegative) integer, and 
f(d 2 ) remains unspecified. Interestingly, the equation becomes much simpler in certain cases. For example, if we 
restrict n to be 0 or 1 , then the last two terms on the right-hand side automatically vanish: when n = 0 the operator 
d x v — vd x vanishes, and when n = 1 it is simply v' and the two terms cancel out. We shall thus restrict ourselves to 
these two cases to exploit this simplification. (We note in passing that the last two terms on the right-hand side of 
(A.3) only contain dervatives of v of at least second order; therefore, in de Sitter space where v is proportional to x , 
they vanish identically for all values of n.) 

It is also useful to choose the parameters of Eq. (A.3) such that the differential operator on the right-hand side 
factorizes (thus leading to the u-v decoupling described in §IIC). This is done on the first line, splitting the d 2 into 
two occurrences of d x which are absorbed into the terms in brackets: 


[(d t + d x (v ± 1 ) + (-id x ) n 7 (-id x ) n ) (d t + (v T 1)3* + (-id x )y(-id x ) n ) 

±(—1 ray ay 1 T yiyayya: + m 2 + / ( -a 2 x ) - (ayay 2 ] </, = ya x ) n j(t , X ) (a.4) 


Note the sign ambiguity, which depends on how we choose to place ±9^ in the two brackets of the first line. The idea 
is to make Eq. (A.4) exactly factorizable by choosing parameters such that the second line vanishes. The first two 
terms of the second line give 

± (- 1 )" (ayay 1 - ayyay = ±(-i ya: ( 7 a x - a xl ) a n x = T (-i ) n ay'a n x . (a.5) 

This vanishes when 7 ' = 0, in which case the two possible orderings of (v + 1) and (v — 1) in the first line become 
equivalent. On the other hand, when 7 ' ^ 0, it is still possible to factorize the equation by choosing m 2 and/or f(P 2 ) 
appropriately so as to cancel out the second line of (A.4); but since this choice depends on the sign of the term in 
(A.5), the ordering of (v + 1) and (v — 1) will matter. 

The simplest way to factorize Eq. (A.4) is to take n = 0 or 1, and, if n = 1, take 7 to be constant. Then, to cancel 
out the second line of (A.4), for n = 0 one should take the bare mass function to be m 2 = 7 2 ± 7 ' and the dispersive 
term / 2 = 0, while for n= 1 one should take ?n 2 = 0 and f(P 2 ) = 7 2 P 4 . (Note that the dimensions of 7 are different 
for the two values of n; in fact, when n = 0, 7 is simply the dissipation rate and is usually called T, as in App. F.) 
Finally, in both cases, we have 


d t + 7 (-idx) 2n +d x (v± 1 ) 


d t + 7 ( -id x ) 2n + (v =F 1 ) d x 


<i> = (-d x ) n J(t,x ), 


(A. 6 ) 


which, with a choice of sign, is exactly Eq. (13). 


Appendix B: Retarded Green function 


Here we generalize the considerations of §111 A, finding the form of the Green function for a general wave equation 
with polynomial dispersion relation, for a system which is asymptotically homogeneous. This can be considered as a 
generalization of the Jost construction [39]. 

Firstly, we consider the purely homogeneous case. Then the Fourier transformed (in both space and time) Green 
function is found straightforwardly from the Fourier transformed wave equation, and is of the form 

G “ J "" (t) = ’ (B,1) 

where the superscript R (L) corresponds to a positive (negative) imaginary part of k. so that the mode decays to the 
right (left). Inverting the spatial Fourier transform, and using the Cauchy residue theorem, we then have 


G:,ho m (v') =iO(x-x')Y J *r*e ik ^-*' ) +i6(x'-x)Y / K eikf{x ~ X ' ) 


(B.2) 
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where we have defined 


M? = 


Go Yii^k ( k k - fc f) n, ( k k - k . 


Mt 


1 

Go n, i k k - k f) u^k ( k k - k f)' 


(B.3) 


The inverse coefficients 1/N R ^ L can be thought of as generalized squared “norms” of the unit plane waves. Note 
that they can be written as derivatives of 1/G“ t hom , which itself is just the dispersion relation. In the relativistic 
and dispersive cases, this means that they can be written as ±2 (ui — vk)v g (k), where v g (k) = cLu/dk is the group 
velocity, see [10]. 

Now consider an inhomogeneous system which becomes homogeneous asymptotically. In the diagonal asymptotic 
regions x, x' —> ±oo, the Green function will be G“ t hom (:r, x') of Eq. (B.2) plus a solution of the source-free form of 
the wave equation, i.e. a solution which is just a sum of products of plane waves containing no Heaviside step functions. 
This additional term must obey appropriate boundary conditions; in particular, it must vanish asymptotically, in both 
x and x '. So, in the right-hand asymptotic region x,x' —> oo, we can add to (B.2) a term of the form 


*EE 


Aky'x — ik 1 ; 

Cije z e 3 


since this is the only source-free solution that converges as x,x' —> oo. Therefore, the general form of the Green 
function in the right-hand asymptotic region is 


G?et(x,x’) = i0(x - x') ' 


ik* 


Mi 


■R £ -ikf 


-iki 




+ i6( x' — x) 


E 




ik • 




—ik 1 ; 


i |_ 3 J 3 L i 

(B - 4) 

Let us emphasise that Eq. (B.4) applies only in the right-hand asymptotic region where the system becomes homo- 

R/L 

geneous. In particular, the wave vectors k k ' are solutions of the dispersion relation in the right-hand asymptotic 

region, and the “normalization” factors M k ^ L come from the expansion of G“ t hom (£, x') using the background values 
in this region. Similarly, the coefficients Cij will be particular to the right-hand side, and are fixed by the require¬ 
ment that, when this asymptotic solution is extended throughout the space, the resulting Gf et (x,x') must vanish 
asymptotically. The functions of x and x' appearing in Eq. (B.4) - which are solutions of the wave equation and 
its dual, respectively - can then be used as a set of basis modes. However, we make another choice here: since the 
interpretation of the anticommutator in §IV A relies on the modes (x) containing a single outgoing wave in one of 
the asymptotic regions, we impose that this be the case. Those modes multiplying 9(x — x') in the first term of (B.4) 
are exactly of this form, and so we define (j\' R (x ) and <j>J’ R (x') such that, as x,x' —t oo, 


±<jj,R / \ 

<Pi (x)^e 


ik. 


7uj , R ( / \ 

0 .’ \x ) y e 


— ik t 


E' 


-ifcf 


'i,3 ' 


(B.5) 


Thus we see that the definition of x) is straightforward while that of (j\’ R {x') is more subtle: it contains a single 
growing mode in the right-hand asymptotic region, but must be purely decaying when continued into the left-hand 
region, and thus contains whatever combination of decaying modes in the right-hand region makes this so. 

The second collection of wave functions in (B.4), which multiply 9{x' — x ), are not of the required form, but can 
be made so by redoing the above analysis in the left-hand asymptotic region. Using bars to denote wave vectors and 
“normalization” factors evaluated on the left-hand side, the homogeneous form of the Green function there is 


G r V om 0r,:r') = i0(x - a/) + i9{x' - 


(B.6) 


to which we can add a term of the form 


*EE' 


Akt 


— ik- 




since this is the only source-free solution that decays as x,x' —> —oo. So the asymptotic form of the Green function 
as x, x' —oo is 


G^et (X,x')=i9(x-X')Y^ 

N R e ik * x + ^ x 

e- ih ? x ' + i0(x' ~x)J2 

M R e~ ih i x ' + ^2 c id e~ ik ? x ' 

i 

3 

3 

i 


(B.7) 
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The second term is now of the form required, and we define <l>. 


L {x) and (fh’ L (x') such that, as x,x' —> —oo, 


±uj, L 

<t>3 


(x) e 


ik^x 


4> U i' L {x') -t e 


—ik^x 


+ 


1 

Mf 


E< 


—iky'x 


(B.8) 


where, as before, the coefficients Cij are fixed by the requirement that <(h ,L (x') be purely decaying when continued 
into the right-hand asymptotic region. 

Finally, using the modes defined in (B.5) and (B.8), the full Green function can be written 


GT e = id{x - x')^,^ 4 >T R {x) 4>i u,R {x') + id{x' - x)^^ 4>j' L (x)^j “ ,L ( x '). 


(B.9) 


Appendix C: Modes and dual modes for a simple black hole profile 

In this appendix, we give explicit expressions for the flow velocity profile and the modes and dual modes used to 
construct the figures in Secs. Ill and IV. 

The adimensionalized velocity profile described by Eq. (30) is taken to be 

h(X) = tanh(X) . (C.l) 

This satisfies the conditions given just after Eq. (30): h( 0) = 0 and h'{ 0) = 1, and it asymptotes to constant limiting 
values in the right (subsonic) and left (supersonic) regions, with flow velocities — 1 + D and — 1 — D, respectively. We 
note that the profile is symmetrical in the sense that \v + 1| is symmetrical around the horizon. 


1. Modes 


Let us consider stationary modes of the form <j)(T,X) = (j> w {X)e ™ T , where w = uj/k is the adimensionalized 
Killing frequency. Then Eq. (34a) becomes 

[—iiv - X 2 d 2 x + tanh (X) d x ] <TP0 = 0. (C.2) 

First of all, we note that since the velocity profile is asymptotically constant, the general solution in the asymptotic 
regions decomposes into a sum of exponentials, and Eq. (C.2) reduces there to a dispersion relation: 

A 2 K 2 + ihK — iw = 0, (C.3) 


where h = 1 in the right-hand (subsonic) region and h = — 1 in the left-hand (supersonic) region. In the geometrical 
optics approximation, we can think of the roots as varying continuously with X as h varies between —1 and 1, and 
it is precisely this migration of the roots that is shown in Fig. 2. Note that what is plotted there is the root with 
positive real and imaginary parts, i.e. it is propagating to the right in the subsonic region (the other root, which 
propagates to the left in the supersonic region, follows curves which are simply the negatives of those in Fig. 2, due 
to the invariance of Eq. (C.3) under h —h and K —> —K). It is precisely the dissipative version of the standard 
right-propagating root; for this reason, we label it while the other root, which decays to the left, is K u l- As 

h (and X) are decreased, the real and imaginary parts of K u r both increase, becoming equal at the horizon where 
h vanishes, and where, according to Eq. (C.3), K oc y/i. As h is decreased further, the real part of K u r decreases 
again while the imaginary part continues to increase, until at h = —1 (X —oo) it settles on a value with the same 
real part it started with but a much larger imaginary part. It has thus, in the left-hand region, become the additional 
root due to the dissipative fc 4 term in the dispersion relation. 

Equation (C.2) is exactly solvable using trigonometric and hypergeometric functions. Firstly, <f> w (X) is written in 
the form 


r(x) 


ijj w (X) exp 



r x 


tanh(X , )cLV , 


ip w {X) (2 cosh(X)) 1/2A2 , 


(C.4) 


with ip w (X) obeying the transformed differential equation 


.w 1 1 

“ 4A 4 + 2W 



sech 2 (X) 


ip w {X) = 0. 


(C.5) 
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The solution of (C.5) is given in §23, Problem 4 of [40]. One of the two final independent solutions is 

Kr{X) = (2 cosh (X)) iKuR 2 F 1 ^-iK uR ,l+iK uL ;l- l -(K uR ~K uL )-^{l-t a nh(X))^ , (C. 6 ) 

where, following the notation of App. B, the absence of overbars indicates that the wave vectors are to be evaluated 
in the right-hand asymptotic region, i.e. with h = 1 in Eq. (C.3). To see that (C. 6 ) obeys the correct boundary 
conditions, recall that the modes (j) u are defined as containing a single outgoing (hence decaying) wave of unit 
amplitude in one of the asymptotic regions. Taking the limit X —» oo, the hypergeometric function approaches 1, so 
4>u R (X) —> exp (i K uR X) as required. The other independent solution is obtained via reflection, X —> —X : 

<% L (X) = 4% R (-x), (c.7) 

which, in the left-hand asymptotic region, is equal to exp (— i K uR X) = exp (i K uR X) , as required (the overbar on 
K u l indicating that it is to be evaluated in the left-hand asymptotic region, with h = — 1 in Eq. (C.3)). 


2. Dual modes 


In constructing the retarded Green function for the <fi field, we need its dependence on both the image coordinate 
x and the source coordinate x '. As explained in §IIIA, the latter obeys the dual of Eq. (13); at the level of the u- 
and u-rnode source-free equations, this amounts to switching the sign of 7 (and as such describes dissipation acting 
backwards in time). Using a tilde to denote the dual mode, it obeys 

[iw + X 2 d 2 x + tanh (X) d x \ 4>~ W {X) = 0. (C. 8 ) 

At the level of the dispersion relation (C.3), this effectively switches the sign of h so that the wave vector solutions 
for the dual wave equation are simply the negatives of those of the original wave equation. This means that the 
propagating solutions now steadily increase in amplitude towards infinity, whereas the additional complex modes 
decay rapidly to zero. 

Equation (C. 8 ) can be solved exactly in a similar manner to Eq. (C.2). One of the independent dual modes is 
0~™(X) = A~ w ( A) (2 cosh (X))~ iKwL 2 F 1 (iK uL , 1 - iK uR -1 - l - (K uR - K uL ) ; i (1 + tanh (X))^ , (C.9) 

where the amplitude A _t "(A) is given by 

A~ w (\\ — _ r ( iK uR ) T(1 - iK uR ) _ 

“ r(i-| (. k uR - k uL )) t (-| (k uR - k uL )) ' ’ 

The amplitude is fixed by the boundary conditions imposed on (j )~ R , which (as mentioned in §111 A) are that it decays 
to the left and has a growing mode exp (—i K uR X) with unit amplitude on the right. The decay to the left means 
that it must be proportional to exp(—iK uR X) = exp (iK uR X) in the limit X —> — 00 ; examination of Eq. (C.9) 
shows that this is the case. To ensure that the growing mode on the right has unit amplitude, we use the following 
transformation law in order to decompose the dual mode into its exponential components on the right-hand side: 


2 F 1 (a,b;c-,z) = — ^ 2 A) (a, b ; a + b + 1 — c; 1 — z) 


T(c — a)T(c — b)' 


T(c)T(a + b — c) 

r(o)T(6) 


( 1 -*) 


z—a—b 


F 1 (c— a, c— 6;l + c — a — 6; 1 — z). (C.ll) 


The amplitude (C.10) thus ensures that the coefficient of exp (—i K uR X) is 1. 

As before, the other independent solution can be obtained via reflection about the horizon, X —> — X : 


CW = GK)- 


(C.12) 


The symmetry properties relating the various wave vectors ensure that <PZl(X) also obeys appropriate boundary 
conditions. 

Examples of the dual inodes can be seen in Fig. 3, since G™ t (X, X') is proportional to them at fixed image 
coordinate X. In the top row, we can see the evolution of a propagating wave into a strongly damped wave as the 
horizon is crossed, and that this gradual transition becomes a sharp cutoff in the relativistic limit A 2 —» 0. 
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Appendix D: Noise kernel 

In this appendix, we derive the form of the noise kernel in Eq. (44) from the action of the free environment field 
given in Eq. (7). 

The noise kernel N UI (x, x') represents the sourcing of the field due to the homogeneous part 4/° of the environment 
field. We thus need an explicit form of the operator 4y) in order to find it. For this purpose, it is convenient to work 
in a coordinate system (r, z) adapted to freely-falling observers, in which the oscillators of the environment are at 
rest: 


r = t, 


dx' 

v(x') 


= t. 


(D.l) 


In this coordinate system, the partial derivative d T = dt + v d x , which is exactly the differential operator d T introduced 
in Eq. (8) and the unit timelike vector field u of Eqs. (2) . The free action Sq, of the environment takes the form 


Sy = — f drdz j dq{ ( d T 


'$*• - b 


bh 


(D.2) 


and the free equation of motion for 4/ becomes 




ll M^o = 0 

v(z) q 


(D.3) 


Restricting 4^ to be real, the general solution is 


4 / o( r > z ) = 


K z ) 

){x) 


i( z ) 


\/2oJ rj 


+ C*Jz)—=) , 


V 2 . 


Ida 


(D.4) 


where the factors of 1/ co q are included to give normalized plane waves of unit Wronskian. From the action (D.2), 
the canonical momentum is 


n u 


v(z) q '' 


K z ) 


(D.5) 


and upon quantization and the imposition of the equal time canonical commutator 

z), n °/( T > z> ) = iS(q- q') 8(z- z') 

we find that the quantum amplitudes c q (z) obey the bosonic commutation relations 


(D.6) 


c q (z),ct(z') 1 = S(q - q') S(z - z'). 


(D.7) 


Considering first the case n = 0, the source operator j is 

J = 7(z) 


V ( Z )_g l V ( X ) 


= —i r ){x) 


v(x) y v(z) 

v{z) 


dq'Z" 


){x) 


I g L '~ u q i g 4 ^ u q 

dqu) q Cg(z) — - cUz) 




Using Eq. (D.7) and the fluctuation-dissipation theorem in a heat bath at temperature Tip, 

c q (z),cl{z')}) =cothf^- N ) \c q (z),cl(z f ) 


= coth 


2 T®, 

(it) 


S(q - q')S(z - z') . 


(D.8) 


(D.9) 
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FIG. 8. Variation with D: Here is plotted the relative difference in n u (solid curves) and n_^ (dotted curves) as a function 
D , the reference values being taken at D = 0.5 (i.e. the values plotted in Fig. 6 ). The temperature is fixed at 7 V/k = 10 ~ 2 
and A 2 = lO ^ 1 (so 7 varies as D 2 ). The left panel has ui/k = 1/2-7T, the lowest of the frequencies plotted in Fig. 6 , whereas the 
right panel has uj/k= 1/2, the highest frequency plotted in Fig. 6 . 


we find 


N(t,x;t',x') = (j J(f,x), J(t , ,a; / )|) 


= 5 \ t — t' — 


dx" 


7(x)7(; 


Sf/^o-h (^) c« (u,l ^y) . (D. 10 ) 


v(x")J s /v{x)v{x') 

Because of stationarity, this is seen to be a function only of the time difference t — t! , and we can Fourier transform 
in time to get 

/ +00 

d(t — ^iV(t, x; t', x') 

-OO 

r x dx " 


7(x)7(x') 


=exp 


r dx" 

x' v ( x ") 


UJ n 


dqu} q coth ( ) cos ( u q 


sjv{x)v{x') 

hak mkh' u ’ Mb (S:) exp ( i( “ - L 


x' v(x") 


2tt Ak /x(x)x(x / ) 

where in the last line we have effected the change of variable u q = irq. 


dx" 

v(x") 


(D.ll) 


Appendix E: Subdominant effects 

In the body of the text, we saw that the main parameter governing the deviations from the non-dissipative case 
is A 2 of Eq. (35). In this appendix we study the residual effects of changing D and increasing the environment 
temperature. We also study the strength of the correlations between u- and x-modes, which are here very small since 
they are induced by the coupling to the environment. 


1. Changing the extension of the near-horizon region 

To examine the effects of D , in Figure 8 we plot the relative change in n u and from their values at D = 0.5, 
at fixed low temperature T^/k = 10 -2 and dissipative parameter A 2 = 10 _1 . There is clearly an increased sensitivity 
to D at higher frequency. Also very clear from these plots is the asymmetry between n^ and n_ w , whose values tend 
in different directions. However, the largest total variation seen in any of the observables is around 10%, while D 
varies by almost 100% (i.e. from 0.05 to 0.95). We can thus conclude that it is the combined parameter A 2 = 7 k/D 2 
which is most relevant in determining the values of the observables, while variation of the individual factors (keeping 
A 2 constant) produces only a subdominant variation. 





















23 



FIG. 9. Observables at high temperature: Here are plotted the same observables as in Figs. 6 and 7, except that the 
temperature is now fixed at the higher value of Tq, / k = ir/2 (about 10 times the Hawking temperature). D is again fixed at 0.5, 
and the measurements are made at fixed distance from the horizon, k \x\ /D = 5. The first row shows the relative difference of 
the subsonic quasiparticle number from the standard Hawking values (note that, in the upper left plot, the vertical axis uses a 
logarithmic scale); again, the corresponding value in the supersonic region behaves similarly, and is not shown here. The dots 
in the upper left plot show the data points of the solid curve corrected for the mode propagation from the emission region to 
the detector. The second row plots the rescaled nonseparability parameter, with the shaded area corresponding to A w < 0 and 
hence to a nonseparable state. In the left column, the variable is w/k, and the various curves correspond to different values 
of A 2 : 0.1 (solid curve), 0.2 (dashed curve), 0.3 (dot-dashed curve) and 0.4 (dotted curve). In the right column, the variable 
is instead A 2 , while the various curves correspond to different values of ui/k: 1/20-7T (solid curve), 1/27T (dashed curve), 1/3 
(dot-dashed curve) and 1/2 (dotted curve). 


2. High temperature effects 


As it propagates outwards from the horizon, the <f> field will thermalize, and if the environment temperature TV 
is sufficiently high this will cause the particle number to approach an equilibrium value n eq . Since the environment 
temperature applies in the freely-falling frame, n eq is determined by the “preferred” frequency f l: 


2n eq T 1 


coth 


n 

27V 


coth 


W - Vas k 

27V 


(E.l) 


where v as is the asymptotic value of v. Hence, in an asymptotic region, the total value of n as a function of position 
is thus approximately given by 

n(x) ~ n eq + 5ne~ 2Ta - s ^ , (E.2) 

where Sn itself is independent of x. (Note that Sn may be negative, in which case n grows with x.) 

In Figure 9, the same observables as in Figs. 6 and 7 are plotted, but the environment temperature has been 
fixed at the higher value of TV/ k = 7t/2 (about 10 times the Hawking temperature). We now see a change in the 
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FIG. 10. u-v CORRELATIONS: Here are plotted the magnitudes of the two u-v correlations, dZ (left panel) and cZ (right panel). 
The geometry is fixed at D = 0.5, and the temperature at TW/k = 10” 2 . The various curves correspond to different frequencies: 
lu/k = 1/207T (solid), 1/2-7T (dashed), 1/3 (dot-dashed) and 1/2 (dotted). These are much smaller than the correlations between 
the two it-modes, which is of order as can be deduced from the smallness of Au/riu in Figures 7 and 9. 


behavior of the spectrum, with the relative difference now increasing both with A 2 and with w, and the low-frequency 
temperature being higher than the standard Hawking value. As before, we take our measurement of the particle 
number at kx/D = 5, which is in the flat asymptotic region, but close enough to the horizon so that the coupling 
to the environment has the least amount of time to act as the quasiparticles propagate outwards from the emission 

region. We also perform the same correction described in §IV C, multiplying Sn by exp ^ dx' 3fc^(a/)^ to account 

for the field-environment interaction that occurs between the emission region and the detector. We adopt the same 
rule for x em ■ as in §IV C, taking it as the point where the integrand of Eq. (46) and Fig. 4 drops to 1/5 of its maximum 
value. The corrected plot (shown as dots) shows that, for small values of A 2 , the increase in n su b at high frequencies is 
largely due to injection of quasiparticles by the environment during propagation, rather than any increase in emission 
from the black hole. From the right upper plot, when sending A 2 —> 0, we see that for all ui/k, the limiting value 
of the mean occupation number in our dissipative model is that found by Hawking. In other words, even if the 
environment temperature Tip is much higher than the Hawking temperature, if the dissipative scale is sufficiently high 
(with respect to that temperature) the outgoing spectrum is still characterized to very good approximation by the 
standard Hawking temperature, as if one were in the Unruh vacuum. 

The trends in the nonseparability parameter A remain similar to the low temperature case: it increases with both 
A 2 and w, though we are now able to see a loss of nonseparability where A = 0. However, the main result of Fig. 9 is 
that there exists, even at this high temperature, a wide window over which A < 0 and where the Hawking radiation 
remains quantum in character. The origin of this is to be found in the fact that, for parameters in this window, there 
is a large redsliift in the temperature in the near-horizon region which takes place before the pairs separate. As a 
result, the state is close to the standard Unruh vacuum when the pair creation takes place. 


3. The strength of u — v correlations 


Another subdominant effect concerns the strength of the correlations between outgoing u- and u-modes. It is 
shown in Figure 10. There are two such correlations, d v u and c“, corresponding to the amplitudes of the terms 
<f% R (x) and (4>~^{x))* (</>“(a/))*, respectively, in the full expression for G“ c (x, x'), see Eq. (38). The strength 

of these correlations is always much smaller than that of the u — u correlations encoded in c w . It is seen to increase 
approximately linearly in A 2 . They behave in opposite ways regarding their w-dependence: as lo increases, |d"| 
decreases while |c“ | increases. It should be emphasized that the strength of the u — v correlations has a severe impact 
on the nonseparability of the u — u pairs, see [26]. We strongly suspect that the robustness of the nonseparability 
obtained in the present model is to a large extent due to the smallness of the u — v correlations. 
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Appendix F: The case T independent of P 

In this appendix, we consider the case where n = 0 in Eq. (13), so that the dissipative rate F is independent of the 
momentum P, and is in fact simply equal to g 2 / 2. Then, after Fourier transforming in time, the mode equation (18) 
becomes 


□£<T = l-i (w + *r) + d x (v + 1)] H (W + iT) + (v-l) d x ] 0“ = 0, (F.l) 

while its dual is 

□ = [i (w + iT) + d x {v- 1)] [i (w + ir) + (v + 1) 8 X \ 4r u = 0. (F.2) 

Note that we can remove the first factor in square brackets from each of Eqs. (F.l) and (F.2) to get simpler, first order 
differential equations whose solutions are automatically also solutions of (F.l) and (F.2). Given the sign appearing in 
the v± 1 term, these solutions are u-modes for (F.l) and dual u-modes for (F.2). As mentioned in §IIC, the ordering 
of v + 1 and v — 1 is irrelevant when T is independent of x, in which case all solutions are solutions of simpler first 
order equations. However, in the general case where T is ^-dependent, this is only true of the solutions singled out 
above. So, while we cannot give an explicit general expression for the it-mode <j>^, we can for its dual: 


where we have defined 


<t> u u { x ) = exp 




w + *r 

V + 1 



exp (— iujx u ) x A. 


x u 



A u = exp 



(F.3) 


(F.4) 


This explains the particular ordering chosen in Eqs. (F.l) and (F.2): since it is the dual mode <(>^ w (x) that is integrated 
over to find the anticommutator, we wish for it to have the simple expression (F.3). Although the w-mode equation is 
first order, it does have two independent solutions due to the divergence of the integrand in Eq. (F.3) at the horizon 
where v + 1 = 0; since F > 0, (j)-^ must vanish there. So, as for the T oc P 2 model studied in the main body of this 
paper, there are two distinct u -modes localized on opposite sides of the horizon, which we label cj) l ^ R and and 
their dual modes <j)~ R and 4>Zl- 

The analysis of the main paper progresses in much the same way, and for the anticommutator, Eq. (46) becomes 


\i (I))) .(X. ’G 1 ;;,(x,xi) ■ y' C ,; (x',x 2 ) j x 7rTi.coth (jiT* (r 0 (xi) - t 0 (^2))) . 


There are two differences with respect to Eq. (46): there are now no derivatives acting on the Green functions, since 
there are no derivatives in the n = 0 version of Eq. (43); and 7 has been replaced by ^/r(xi)r(x 2 ) = g(x\)g(x2)/2, 
since g is not now assumed to be x-independent. The reason we have allowed g , or equivalently F, to depend on 
x can be appreciated by examining Eq. (F.5) and comparing it with Eq. (46). We expect the main contribution 
to the integrand for G“ c to come from the near-horizon region. In the F oc P 2 case, Eq. (46), this is achieved via 
the appearance of the x-derivative of Gf et , which is largest around x = 0. In Eq. (F.5), however, setting T to be 
constant everywhere means that all regions contribute equally, and it becomes difficult to distinguish the contribution 
from the black hole and that from the local environment. Another way of saying this is that, while the corrections 
to the Hawking spectrum should be of order O (Thor/w), where Thor is the value of T in the near-lrorizon region, the 
quasiparticle number in the asymptotic region is only well-defined up to a correction of order O (r as /w), where r as is 
the asymptotic value of T; therefore, if T^or = r as , it is not possible to give a meaningful correction to the Hawking 
spectrum due to dissipative effects. We thus introduce an asymptotic switch-off of T, so that r as = 0; yet, while 
this means that the asymptotic fluxes are unambiguously defined for any given v and T profiles, it also introduces 
extra degrees of freedom in the precise form of F. We demand that T(x) satisfy two criteria: it is constant in the 
near-horizon region, and it falls off to zero at a rate slower than k, so that spectral corrections due to a non-adiabatic 
switch-off are negligible. 
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R Q [2oj G“ Uu (x,x')] 



FIG. 11. Green function with P-independent T: Here is plotted, for F independent of the momentum, the w-part of 
the retarded Green function, G“ et!U (ai, x') as a function of x' with x fixed (solid curve). For convenience we have used the 
transformed (tortoise) coordinate x u of Eqs. (F.4), which in the near-horizon region behaves like x u ~ log( ie)/k. We have 
taken uj/k = 1/27T, D = 0.2, and x such that x u = 100. The dashed line plots T/Thor, where Thor is the value of V at x = 0; here, 
we have taken rhor/K = 10“ 2 . The dotted lines plot dcA u of Eqs. (F.4), and since this is w-independent, changing the frequency 
would change the wavelength in x u but not the envelope, which behaves like x rhoi ^ K = exp (Fhor*u) in the near-horizon region. 


Results for a low environment temperature T\j,/k = 10 -2 are shown in Figure 12. Interestingly, just as for the 
r oc P 2 model, we find that the corrections to the Hawking spectrum amount to a global change of temperature in 
the low-frequency regime u>/n <C 1, where the relative change in the spectrum (in fact a reduction) is proportional 
to ri lor /«;. The coefficient is approximately 20. Thus the robustness of the spectrum is again recovered when the 
following inequality is satisfied: Thor/ft < 1/20. 

The differences from the T oc P 2 model show up clearly when we look at the impact of dissipation on the strength 
of the correlations. When studying the parameter A u , we find that the pairs of tt-quanta are nonseparable at high 
frequencies and separable at low frequencies, with the changeover occurring around u> ~ Th or - This result is relevant 
for experiments in polariton systems which aim to observe the spontaneous channel [17, 18]. 
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